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Grid Generation About Complex Three-Dimensional 
Aircraft Configurations 

The study of the aerodynamics or radar cross-section signatures about complete 
aircraft by numerical methods is still limited by two (or more) problems. The 
first involves the lack of adequate turbulence models applicable to a wide range of 
flow regimes and the second has to do with obtaining three dimensional grids with 
sufficient resolution to resolve all the flow or other physical features of interest. This 
study is concerned with the latter problem. 

The generation of a computational grid involves a series of compromises to 
resolve several conflicting requirements. On one hand one would like the grid to be 
fine enough and not too skewed to reduce the numerical errors and to adequately 
resolve the pertinent physical features of the flow field about the aircraft. On the 
other hand the capabilities of present or even future supercomputers are finite and 
the number of mesh points must be limited to a reasonable number, one which is 
usually much less than desired for numerical accuracy. 

One technique to overcome this limitation is the ’zonal’ grid approach. In 
this method the overall field is subdivided into smaller zones or blocks in each of 
which an independent grid is generated with enough grid density to resolve the flow 
features in that zone. The zonal boundaries or interfaces require special boundary 
conditions such that the conservation properties of the governing equations are 
observed. Much work has been done in 3-D zonal approaches with nonconservative 
zonal interfaces. A 3-D zonal conservative interfacing method that is efficient and 
easy to implement was developed during the past year. 

As originally envisioned, the conservative interfacing was to be developed for 
finite difference codes which were the most prevalent codes at NASA/Ames. How- 
ever during the course of the work it became apparent that such a procedure would 
become rather complicated and that it would be much more feasible to do the con- 
servative interfacing with cell-centered finite volume codes. Accordingly, the CNS 
code was converted to finite volume form. This new version of the code is named 
CNSFV. The original multi-zonal interfacing capability of the CNS code was en- 
hanced by generalizing the procedure to allow for completely arbitrarily shaped 
zones with no mesh continuity between the zones. While this zoning capability 
works well for most flow situations it is, however, still nonconservative. The conser- 
vative interface algorithm has also been implemented, but has not been completely 
validated yet. 

The results of this work has been reported in reference 1, a paper to be pre- 
sented at the 10th AIAA CFD conference to be held 24-27 June 1991. This paper 
essentially covers all the details of the work done under this contract and is included 
as an appendix. 
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Conservative Multizonal Interface Algorithm 
for the 3-D Navier-Stokes Equations 
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Abstract 

One method of solving the Navier-Stokes equations 
about complex and realistic aerodynamic configurations is 
to use a zonal method. In this method the overall flow 
field domain is subdivided into smaller blocks or zones. In 
each of these zones, the flow field is solved separately of 
the other zones. The boundary data for each zone is pro- 
vided by the neighboring zones. The major difficulty of 
the zonal methods has been how to maintain overall con- 
servation for arbitrarily shaped zones A new method of 
conservative patched zones has been developed. It uses 
structured meshes in the individual zones. The interface 
between the zonal block faces is defined by the union of 
the face points of adjoining blocks. An unstructured gTid 
is generated upon which the interface fluxes can be e- 
termined. Flux balancing of the interface fluxes is then 
easily achieved to obtain global conservation. The method 
has been implemented into two Navier-Stokes codes. The 
use of the procedure is easily implemented into other fi- 
nite volume codes. There are no topological restrictions 
on the zonal boundaries; e.g. the zonal interfaces can be 
curved surfaces for ease in constructing structured meshes 
in each of the zones. Several examples axe presented to 
demonstrate the viability of the interfacing procedure. 


Introduction 

There are two basic approaches of numerically simu- 
lating the Navier-Stokes equations about complex and re- 
alistic aerodynamic configurations. One is based on struc- 
tured meshes in which the neighbors of a mesh point are 
known implicitly. The other approach is based on an un- 
structured mesh in which the neighbors of a mesh point are 
not known implicitly and this information must be stored 
for each point. Numerical methods based on structured 
meshes are well developed, but suffer limitations when 
dealing with complex configurations in that it is difficult, 
if not impossible, to generate a single mesh about such a 
configuration and still have the required mesh qualities for 
stable and accurate numerical solutions. The generation 
of unstructured grids about complex configurations is, in 
principle, much easier; however, the numerical methods for 
such grids are not yet mature enough to compete against 
structured grid numerical methods. 
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The structured grid generation problem for complex 
configurations can be alleviated by zonal methods. In 
this method the overall domain is subdivided into smaller 
blocks or zones. In each of these zones, a grid is gener- 
ated and the flow field is solved independently of the other 
zones. The boundary data for each zone is provided by 
the neighboring zones. The major difficulty of the zonal 
methods has been how to maintain overall conservation for 
arbitrarily shaped zones. The resolution of this difficulty 
is the purpose of this paper. 

There are two types of zonal methods in common us- 
age today. One is the overlaid zones in which both zones 
share a common interface region, as used in the Chimera 
[1] approach. The other is the patched zonal technique 
where the two zones share only a common boundary. Ex- 
amples of this approach are the zonal method of Rai [2] and 
Thomas [3]. Both methods have advantages and disadvan- 
tages depending on the application. With patched gnds it 
is much easier to maintain conservation since there is o y 
one boundary across which the two zones communicate. 
On the other hand for moving multiple bodies the overlaid 
grid approach has certain advantages if conservation is not 
important, Le., the flow field is continuous with no shock 
waves or shear surfaces. If conservation is important, as it 
is in this investigation, then we are restricted to patched 
zones. 

Most methods using conservative patching methods 
are restricted to interface surfaces which are planar sur- 
faces due to the minor gaps and overlaps that occur at a 
curved interface if the two zones are not mesh continuous. 
Furukawa et al [4] attempted to resolve this problem by 
using only one of the zones to determine the zonal bound- 
ary for both zones. The open question that remains is 
which zone determines the interface boundary. Furukawa 
et al chose to use the zone which has the better resolution. 
This can result in loss of accuracy if the mesh ratio changes 
in the interface, for example at a viscous boundary layer. 


In this paper the zonal interface boundary is deter- 
mined by the union of all the face mesh points of both ad- 
joining zones. The interface surface is now unique and de- 
termined much more accurately than either one of the indi- 
vidual face surfaces. The collection of interface points is in 
general no longer structured and readily available unstruc- 
tured grid generation techniques can be used to construct 
(triangulate) an interface grid. With the triangulated in- 
terface grid, the metrics (i.e., surface area normals) and 
cell volumes can be determined for each of the interface 
cells. 

The development of the interface algorithm is dis- 
cussed in the following sections. The procedure hai 
been implemented into two finite volume three-dimensional 
Navier-Stokes codes, namely TUFF [5] and a finite volume 
version of CNS [6] henceforth called CNSFV. We will limit 
the exposition to the diagonal version of the Beam and 
Warming scheme [7] as used in the CNSFV code, but re- 
sults from the TUFF code will also be presented. 
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Numerical Scheme 
Navier-Stokcs Equations 

The three-dimensional thin-layer Navier-Stokes equa- 
tions in strong conservation law form in curvilinear coor- 
dinates are 

drVQ + dtF + dvG + d^H =Rc~ 1 d<S (1) 


The metrics used above have a different meaning for a 
finite volume formulation compared to the finite difference 
formulation of '6j. Referring to a typical finite volume ceU 
as shown in figure 1, the finite volume metrics are defined 
(see, for example, Vinokur [8]) as 
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The contravariant velocity components are defined as 
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The finite volume metrics represent the cell face area 
normals in each of the curvilinear coordinates ({, *?.<)• 
They are related to the metrics introduced in equations 
(1 - 5) as follows 


and the viscous flux is given by 
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The pressure is given by the equation of state 
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The volume of the computational cell is given by 


d Q G = b = r„A,r n 1 


V = -[(r 4 - ri) x (r s - ri) • (r T - r s ) 

6 
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and is the finite volume equivalent of the inverse Jacobian 
of the coordinate transformation in the finite difference 
formulation of [6]. 

Diagonal Beam- Warming Algorithm 

The implicit Beam- Warming scheme for the finite 
volume formulation is given by 


[7 + + M n + hC n - Rt~ l S<M n )]SQ- 1 - IT 

where 


R n = _ 4- 6r,G n 4- 5<II n - Re-'ScS") (11) 

The convective three-dimensional flux Jacobians A, 
B, C and the viscous flux Jacobian are defined in the ap- 
pendix of [5]. With the use of approximate factorization 
and diagonalization of the flux Jacobian matrices, a scalar 
pentadiagonal algorithm [7] can be derived as 


T([I 4 V h8^\^]N[I 4 * V 1 h8 n \ 1i }P- 

[I 4 V-^hS^T^ AQ n = R n (12) 

where 6{ is a central difference operator and A<? n = 
Qn+] _ Qn Qn+1 _ + The viscous terms 

are not included in the left-hand (implicit) side. The ar- 
tificial dissipation is included in both sides and is derived 
below. 

The inviscid flux Jacobians are diagonalized as fol- 
lows: 

d Q F = A = T t K t Tf l 


dqE =C = T < \<T ( ~ l 

The Tt,T v ,T( are the eigenvector matrices of A, B, C, 
respectively with A {t A„A ( as the respective eigenvalues. 

We also have 

N = T i ~'T v 

p = r,- , r ( 

Each of the factors of the implicit operator of equa- 
tion (12) has an artificial term added to stabilize the cen- 
tral difference operator. The dissipation is based on Jame- 
son’s nonlinear second and fourth order dissipation and for 
the 4— operator takes the form 

hV~' Vf L*;+i( e<,) ■ - « (4) 01 AQ" (13) 
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where k 2 , k* are constants of o(l), and A{, V( are t l* e 4" 
ward and backward difference operators. & j+ i is a modi- 
fied spectral radius defined as 

*j+i = &j + &j + J 

ffj = cry (l + (17) 

<7j = [|lf|+a v /sx I + s v s + ^ J ]. (18) 

and where cell centered surface areas are used, e.g. 

3xj = ^(•Vr+J +- 5 x.>-i) 

The modified form of the spectral radius, equation 
(17), is suggested by Turkel [9] to account for large aspect 
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ratio computational cells as for example in a viscous layer. 
Similar dissipation terms are obtained for the q- and C~ 
operators. The dissipation terms added to the nght hand 
side of equation (10) are identical to those given above 
except that A Q n is replaced by Q n . 

To show that the above scheme is conservative and 
has the telescoping property, the right hand side of (12) 
can be written as 


jt— 


+£,+ > - - Re ^r+i + Re 
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where F,G,H,S are the numerical fluxes and are defined, 
for example for F , as 
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Similar terms are obtained for the numerical fluxes 
in the other two coordinate directions. 


Bo undary Conditions 

To complete the equation set, boundary conditions 
must be specified. With the use of curvilinear coordinates, 
the physical boundaries have been mapped into compu- 
tational boundaries, which simplifies the application of 
boundary conditions. The boundary conditions to be im- 
plemented for external viscous or inviscid flows include (1) 
inflow or far field, (2) outflow, (3) inviscid and (4) viscous 
impermeable wall, and (5) symmetry conditions. For ex- 
ternal three-dimensional flow fields about closed bodies, 
the topology of the grid usually introduces (6) grid singu- 
larities which require special boundary conditions. The use 
of zonal methods can avoid the generation of grid singu- 
larities, but requires (7) special zonal interface boundary 
conditions. For compressible flows these zonal boundary 
conditions should be conservative to maintain global con- 
servation. 

In the finite volume approach, the specification of 
boundary conditions reduces to specifying the appropriate 
numerical fluxes at the boundaries. The details of im- 
plementing boundary conditions (1) through (5) are well 
known and are given in [5] and [10]. The grid singularity 
boundary condition is described below and the interface 
boundary conditions are given in the next section. 


The grid singularity boundary concision is similar to 
the symmetry boundary condition for the inviscd and vis- 
cous fluxes in that there is no flux through that boundary. 
If however, that is all that is done the results shown in 
figure 2a are obtained. These results are the density con- 
tours of a Mach 8 viscous blunt body flow. As shown m the 
figure, a nonphysical behavior appears at the singular line. 
The nonphysical results are due to a local violation of the 
entropy condition [11]. For central difference schemes, the 
artificial dissipation is the only stabilizing (entropy pro- 
ducing) mechanism available. At the grid singularity, the 
spectral radius, and hence the artificial dissipation, van- 
ishes due to the vanishing of the metrics. The introduction 
of Harten’s entropy correction [12] resolves the difficulty 
and the results are shown in figure 2b. 


Interface Method 

Most methods using conservative patching methods 
are restricted to interface surfaces which are planar due 
to the minor gaps and overlaps that occur at a in ’ 

terface if the two zones are not mesh continuous. This is 
demonstrated in figure 3a. Furukawa et al [4] attempted to 
resolve this problem by using only one of the zones to de- 
termine the zonal boundary for both zones. This is shown 
for the two-dimensional example in figure 3b. The open 
question that remains is which zone determines the inter- 
face boundary. Furukawa et al chose to use the zone with 
the better resolution. This, however, can result in loss of 
accuracy if the mesh ratio changes in the interface as for 
example at a viscous boundary layer. 

In this paper the zonal interface boundary is deter- 
mined by the union of all the face mesh points of both ad- 
joining zones as shown in figure 3c for the two-dimensional 
case. The interface surface is now unique and determined 
much more accurately than by either one of the individ- 
ual face surfaces. This presumes that the individual points 
from both zones lie in the interface surface and that the 
interface surface is itself smooth and continuous. The col- 
lection of interface points is, in general, no longer struc- 
tured, and readily available unstructured grid generation 
techniques can be used to construct (triangulate) an inter- 
face grid. With the triangulated interface grid, the metrics 
(i e., surface area normals) and cell volumes can be deter- 
mined for each of the interface cells, i.e., the computational 
cells that touch the interface surface. 

The redefinition of the interface surface has modified 
all the interface cells and they are no longer hexahedral, 
but rather multifacetted. Determining the surface area 
normals and cell volumes now becomes more complicated. 
There are three types of interface cells, namely, face cells, 
edge cells, and corner cells. The face cells have only one 
face bordering other zones, whereas the edge and corner 
cells have two or three faces in contact with other zones, 
respectively. A typical cell at the interface is shown in 
figure 4. 

The modification introduced by the interface surface 
requires that the metrics and cell volumes of the interface 
cells be corrected to account for the changed shape of the 
interface cells. Corrections are required for the area nor- 
mals of the cell face touching the interface surface, the four 
sidewalls, which may no longer be quadrilaterals, and the 
cell volume. For inviscid steady flows, the cell volume cor- 
rections are not needed since the volume has no effect on 
the steady solution as shown by equation (11). However, 
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for all viscous flows and inviscid unsteady flows the volume 
corrections are necessary. 

There are two ways that the interface grids can be 
triangulated. The first involves eliminating all the grid 
lines from the structured grid cell faces and constructing 
the triangulated mesh. Since in general the new grid lines 
will not be aligned with the original grid lines, a clipping 
algorithm is used to clip those triangles which lie outside a 
particular individual face cell of either zone. An example of 
this triangulation procedure is shown in figure 5. Here two 
zones have square faces each consisting of a Cartesian gnd 
of 15 x 15 uniformly spaced points. The two faces are on- 
ented at an angle of 45° to each other (see figure 5a). The 
union of both sets of face points results in an unstructured 
collection of points which are then triangulated with an 
advancing front unstructured mesh generation procedure 
[13]. The resulting mesh is shown in figure 5b. 

The second procedure retains the original grid lines 
and triangulates any set of points which form a polygon of 
more than three sides. This approach avoids the use of a 
clipping algorithm since each face cell contains an integral 
number of interface triangular cells. An example of this 
procedure is shown in figure 10b for two zones with polar 
grids. 

The unstructured interface grid requires that a set 
of pointers be defined. These pointers indicate which two 
interface cells share a common interface area or section. 
These sections need not be triangular even though the sur- 
face has been triangulated. A common section (or poly- 
gon) may be composed of several triangles if the surface 
is triangulated by the second procedure described above 
or it may be composed of several dipped triangles if the 
first procedure is used. For a more efficient interface flux 
computation it is convenient to define another set cro fs- 
reference pointers. These cross-reference pointers identify 
which interface polygons are in contact with each of the 
interface cells. 


Numerical Interface Flux 

The final step required for the conservative interface 
algorithm is to determine the numerical fluxes at the iuter- 
face. Originally it was anticipated that the interface flux 
could be determined by any stable numerical scheme, not 
necessarily the one used in the interior of the zone. How- 
ever, it was found that if the numerical schemes differ, then 
’’glitches” always appear at the zonal interface. It was also 
necessary to maintain the same order of accuracy for the 
interface scheme as for the interior scheme. An example 
of this case is shown in figure 6. Here the interface flux is 
determined by a first order upwind scheme and the interior 
by a second order upwind scheme. The unfortunate results 
are the discrepancies at the interface boundary as shown 
by the pressure contours. The use of the same second or- 
der upwind scheme for the interface flux eliminates these 
discrepancies. 

The above example shows that it is not possible to 
completely generalize a conservative interface algorithm. 
The numerical flux must be determined by the same nu- 
merical scheme as used in the interior of the zone. The 
geometric aspects of the interface algorithm can be gen- 
eralized, which is essentially the most difficult part of the 
interface procedure. The numerical fluxes at the inter- 
face must (and should) be computed by the particular flow 


solver involved. 

The computation of the interface numerical flux is 
relatively straightforward. Since it is scheme dependent it 
is given is general form first, valid for many schemes up 
to second-order accuracy. A specific form will be given for 
the CNSFV code. Figure 7 shows two zones, jl and j2, and 
the interface with the individual surface polygons labelled 
by V. The flux at the interface section "i" is given by 


F jlJ2 =i = 

where $ is the function defining the particular scheme un- 
der consideration, P is the area normal of the interface 
section w i”, and the ft j and ft 2 are the cell-centered values 
in zones 1 and 2, respectively. 

The numerical flux at the opposite face of the inter- 
face cell is slightly more complicated. It is 


F - 2 _| = jft'2=2»|j2 = 3>ftl = l] 


where ft is the area (U., \P\) weighted average of all 
the interface cells of the zone ”1” sections that comprise 
the face cell of zone ”2” touching the interface surface. 
If higher order schemes are considered, then special care 
must also be taken for the fluxes at the next level of cell 
faces, ie. The above formulation is valid for all 

cell-centered finite volume schemes. 

For the particular case of the CNSFV code, the in- 
terface numerical fluxes arc 


F^ = F^ + a j(e (1) A( Q" - ( (4) A £ V( A ( Q n )j 


where 


F" = i [a'i(/j2=i + /ii=i)” + 3 'v(9j2= i + 9ji=i) n 


and 


+ ^; 2=1 

The flux at j‘2 = | is similarily determined, taking care to 
use the appropriate area averaged values of all cell- centered 
values. 

Results 

To validate the interface algorithm, several tests were 
conducted. The first test is the freestream preserving test. 
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In this test, the inflow and and permeable wall boundary 
conditions (ie. conditions (1).(3) or (4) from the boundary 
conditions section) are set to the freestream conditions and 
the solution is converged. If the initial conditions were also 
set to the freestream condition, then the residual should 
be at machine zero (R= 0(1O“ 14 )) and remain there for 
all subsequent iterations provided that the flow field is dis- 
cretized properly with no gaps or overlaps in any of the 
computational cells and interface boundaries. Indeed, for 
the CNSFV code, this was the case. 

However, the freestream test is not a good indication 
of the accuracy of the scheme. Because of the telescoping 
property of the scheme, the surface area normals can be 
computed inaccurately (even erroneously) and the scheme 
will still pass the freestream test. To test for accuracy, 
a single cell residual, equation (19), is computed with the 
same freestream conditions imposed as above. In this case, 
the maximum residual was of the 0(10 ) on a 64 bit 

machine (Cray YMP). If the grid and the metric terms 
axe computed in double precision (ie. 128 bits), then the 
maximum residual reduces to O(10" 12 ). This indicates 
that machine roundoff errors are not yet a problem, but 
can be if the meshes are refined much further. 

Three different flow cases covering the entire Mach 
regime from incompressible, to supersonic, and to hyper- 
sonic flows with finite rate chemistry were computed with 
the TUFF code. The same basic conservative interface al- 
gorithm described above was used in all three cases, how- 
ever the conservation law equations differed for each of the 
three cases. 

The first case involved an incompressible inviscid flow 
about a cylinder. The two zone mesh is shown in figure 8a. 
The pressure contours are given in figure 8b and the surface 
pressures in 8c. The results across the zonal boundaries are 
smooth and continuous. 

The second set of results are for supersonic blunt 
body flow. This case is an inviscid Mach 2, axisymmetric 
blunt body flow computed on the four-zone mesh shown in 
figure 9a. Figure 9b shows the solution and the bow shock 
position. For these results, the analytic shock location at 
steady state [14] is shown by the solid squares. The solu- 
tion is shown in terms of mach contours on a background 
grid which is cell centered for plotting purposes only. As 
can be seen the multizonal computed and analytic shock 
shapes compare quite well. 

The final flow results obtained are for a viscous hy- 
personic flow about a hemisphere at jV/oo = 15.3 and 
Re = 2.2 x 10 5 /m. The interface triangulation for the 
two zones containing the grid polar singularity is depicted 
in figure 10b. The flow results in terms of Mach and 
atomic oxygen concentration contours are shown in figures 
10c and d, respectively. Again, the results indicate that 
the solution contours are smooth and continuous across 
the zonal boundaries. Although not shown, the computed 
shock stand-off distance agreed well with the experimental 
data of reference 15. 


Closing Remarks 

A conservative zonal interface algorithm has been 
presented. It uses some of the best features of both the 
structured and unstructured mesh CFD technology. The 
interface surface grid is unstructured from which the met- 
rics and interface fluxes can be readily constructed to ob- 


tain the proper conservative interface algorithm. For effi- 
ciency and rapid convergence, the flow solver within each 
of the zones is based on structured mesh CFD technolop. 
The ordering between the zones can be unstructured for 
maximum flexibility in constructing zones and grids about 
complex and arbitrary configurations. The interface algo- 
rithm has been implemented into two three-dimensional 
Navier-Stokes finite volume codes (TUFF and CNSFV) 
and has shown to yield good results. Further testing is 
being conducted for more complex and realistic aerody 
namic configurations. The procedure is general and can 


im 


tr - ' y 

intn nther finite volume codes. 
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a) grid singularity problem. 



b) grid singularity problem resolved 
by Harten’s entropy correction [12]. 


Figure 2. Density contours for a viscous blunt body flow 

M x = 8.0, Reo = 2.19 x 10 5 ,a = 5.0° 
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a) "classical” 

works only for planar interfaces, 
otherwise conservation errors occur. 
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volume correction 



5ide wall metric 
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Figure 4. A typical zonal interface finite volume cell. 
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b) Furukawa et al [4] method 
uses boundary of finer mesh 
zone to determine the interface boundary. 



Figure 5a. The two structured face grids 
of a typical interface surface. 



Figure 3. Zonal Boundary on a Curved Interface 


Figure 5b. The unstructured grid 

of a typical interface surface. 
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Pressure Distribution on a Cylinder 




Pressure contours 


c) Comparison with analytic data. 


Figure 8. Inviscid incompressible flow about a cylinder. 



a) Four zone grid. 


b) Mach contours. 

■ - analytic bow shock location 
(Lyubimov, NASA TT-F715, 1973). 


Figure 9. Axi symmetric inviscid blunt body flow at Moo - 2.0. 
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